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ABSTRACT 

Stochastic heating of small grains is often mentioned as a primary cause of large in- 
frared (IR) fluxes from star-forming galaxies, e.g. at 24 /im. If the mechanism does 
work at a galaxy-wide scale, it should show up at smaller scales as well. We calcu- 
late temperature probability density distributions within a model protostellar core for 
four dust components: large silicate and graphite grains, small graphite grains, and 
polycyclic aromatic hydrocarbon particles. The corresponding spectral energy distri- 
butions are calculated and compared with observations of a representative infrared 
dark cloud core. We show that stochastic heating, induced by the standard interstel- 
lar radiation field, cannot explain high mid-IR emission toward the centre of the core. 
In order to reproduce the observed emission from the core projected centre, in partic- 
ular, at 24 ftm, we need to increase the ambient radiation field by a factor of about 70. 
However, the model with enhanced radiation field predicts even higher intensities at 
the core periphery, giving it a ring-like appearance, that is not observed. We discuss 
possible implications of this finding and also discuss a role of other non-radiative dust 
heating processes. 

Key words: stars: formation - Stars, infrared: ISM - Sources as a function of wave- 
length, radiative transfer - Physical Data and Processes 



1 INTRODUCTION 

Thanks to IRAS, ISO, and, in particular, the Spitzer Space 
Telescope, mid-infrared emission now attracts significant at- 
tention both at small s cale, as an indicat or of the pro- 
tostar formation (e.g., iRagan et al.l [2007), and at large 
scale, as an indicato r of th e global star formation rate (e.g., 
ICalzetti et all 120071 . l2010h - In both cases making reliable 
conclusions is only possible with a radiation transfer (RT) 
model and a detailed account of the dust thermal balance. 

Numerous models have been developed to describe 
emergent spectra for a dusty medium. Most of them treat 
dust radiative heating as a continuous process that leads 
to a common equilibrium temperature for dust grains 
jEgan et al.lll988l: Ihvezic fc Elitzurl Il997l: IWolf et al.lll999l ; 



iDullemond fe Dominikl 120041 : iRobitaille et all 120061 ). How- 
ever, it has been recognized long ago that grains of different 
sizes respond differently to a photon absorption, depending 
on th e ratio of the pho ton energy and the grain thermal en- 
ergy l|Greenberell 19681 ). If the former is greater or compara- 
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ble to the latter, the temperature evolution of a single grain 
is stochastic and consists of short temperature spikes upon 
a photon absorption, followed by prolonged "cold" periods, 
when grain temperature is lower than the equilibriu m tem- 
perat ure it would have if heating were continuous (|Dulevl 
Il973l ). Because of the energy relation, mentioned above, 
stochastic heating depends on the grain size and is only im- 
portant for very small grains (VSG), polycyclic aromatic 
hydrocarbons (PAH), or other similar particles. The over- 
all result for a particular grain is excess emission at shorter 
wavelengths (due to temperature spikes) and lower emission 
at longer wavelengths (due to cold intervals). When a grain 
ensem ble with an MRN-like size distribution (|Mathis et al.l 
Il977l ) is considered, only excess mid-IR emission is seen , as 
emission at longer wavelengths is dominated by large grains 
that are not susceptible to stochastic heating. 

The importance of stochastic heating is well recognized 
by the community dealing with the galaxy-wide star for- 
mation. It now becomes a standard component of galactic 
models to explai n the emergent spectral energy distribu- 
tion (SEP) (e.g.. iDraine fc LillgOOfl ICompiegne et al.ll201~il ; 



tion (StiU) (e. g., Drame & Li 2UU(; Uompicgnc et al.ll2Ulll 
iPopescu et al.ll2011 ; Baes et al. l201ll ). There are also mod- 
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els of protostellar objects and cir cumstellar discs, which in- 
clude the stochastic heating (e.g.. iManske fc Henninell 19981 ; 
iFerland et all 1 19981 ; IWood et al.ll200Sl ). However, in studies 
of individual protostellar objects this process is often ig- 
nored, and 24 /tm emission is rather assumed to be an in- 
dication of the presence of an internal heating sou rce (e.g., 
iBeuther fc Steinacker! [20071 ; iRathborne et alj|20ich . This is 
probably justified in cases when a compact source is seen 
on a 24 (jm image, but the interpretation of diffuse 24 )im 
emission can be less straightforward. 

A common approach to the SED modelling of star- 
form ing regions is an application of the modified black-body 
law (Hil debrandlll983j ). This law is extensively used in the 
analysis of far-infrared, submillimetre, and millimetre obser- 
vations. However, its application is limited only to the emis- 
sion of dust distributed along the line of sight. Accordingly, 
it only allows to find the dust column density and density- 
weighted temperature, which is not enough for the chemical 
modelling needed to interpret line observations. Also, infer- 
ences, based on this law, can be inaccurate and ambiguous. 
For example, observed variations in the dust emissivity index 
P can be caused either by real changes in dust properties, 
by observational noise, or by the temp erature gradient along 
the line of sight (jShettv et al.ll2009al lbl). 

Infrared Dark Clouds (IRDC), that are believed to be 
massive counterparts of low-mass prestellar cores, are simul- 
taneously seen in emission (millimetre and submillimetre) 
and absorption (near-IR and mid-IR). Thus, their proper RT 
modelling hopefully makes the derived density and tempera- 
ture structure less ambiguous and more reliable. Density and 
temperature distributions in a core affect both the shape and 
the depth of the intensity minimum at shorter wavelengths 
and the intensity maximum at longer wavelengths. Thus, to 
reproduce the IRDC shadows and their millimetre emission 
counterparts, one needs a spatially resolved, at least, ID 
model of the RT in the cloud. 

It needs to be taken into account that extended IR emis- 
sion from an infrared dark core consists of three compo- 
nents, namely, attenuated background emission, proper core 
emission, and foreground emission. Proper core emission is 
caused by dust grains heated both by an internal source 
and by the external radiation and can be estimated using 
the core RT model. A major problem of infrared studies is 
the problem of foreground subtraction. There are two possi- 
ble approaches to this problem. The first is to assume that 
the observed intensity at the darkest spot of the considered 
region is entirely caused by foreground emission of any na- 
ture, like zodiacal light, foreground interstellar matter, or 
instrument noise (|Stutz et al.ll2009l ). The second is to take 
the zodiacal l ight contribution fr om the interplanetary dust 
cloud model i|KelsaH et al.lll998h . then to subtract it from 
the signal, and to assume that all the remaining intensity 
comes from the source itself. 

The first approach is safer as in this case one definitely 
knows that derived results represent some limiting values. 
On the other hand, its usage implies an assumption that 
the core is optically thick in the considered range and does 
not produce any proper emission. These are two additional 
and maybe unwanted constraints to the model. The second 
approach seems to be more honest, but it brings up a ques- 
tion of a foreground subtraction. 

The foreground problem can be related to the ex- 



ces s 24 /xm emission, supposed ly found in two IRDC cores 
by IPavlvuchenkov et alj l|201lf) . They attempted to model 
SEDs of infrared dark cloud cores IRDC 320.27+ 029(P2) 
and IRDC 321.73+005(P2) JVasvunina et all I2009T) in the 
range from 8 /im to 1.2 mm. Pavlyuchenkov et al.l (|201lh 
found that a model of a spherically symmetric core, that only 
accounts for equilibrium dust thermal emission, allows to re- 
produce absorption of background radiation at 8 /im, emis- 
sion (or the lack of emission) at 70 fim and 1.2 mm, but fails 
at 24 fim. Absorption minima at this wavelength, predicted 
by the model, turned out to be mu c h dee per than is ac- 
tually observed. IPavlvuchenkov et alj l|201lf ) suggested that 
an excess emission at 24 /im can be generated by stochasti- 
cally heated VSGs, that are w ell-known pote ntial sources of 
emission in the mid-IR range (Kruegel 2003). 

In this paper we use an ext e nded version of the RT 
model from IPavlvuchenkov et all l|201ll ) to simulate emis- 
sion of a dense interstellar clump heated from outside by 
the interstellar UV radiation and possibly from inside by 
a protostellar object. Emission from stochastically heated 
VSGs and PAHs is taken into account. Our initial goal was 
to reproduce shallow 24 / j,m sh adows in IRDC cores studied 
bv lPavlvuchenkov et all (|201lh . While pursuing this goal, we 
found that a stochastic dust heating model of a kind, that is 
used in galactic SED modelling, being applied to individual 
clumps, predicts their distinct morphological features. We 
discuss possible implications of this finding and also con- 
sider a role of other non-radiative dust heating processes. 

The structure of the paper is the following. In Section 2 
we discuss the problem of foreground emission and show 
how its level of uncertainty affects the derived core param- 
eters. The protostellar core model and the RT method with 
stochastic heating algorithm are described in Section 3. Also, 
in this section results of protostellar core simulations with 
stochastically heated dust grains are presented. In Section 4 
we consider possible ways to solve the excess mid-IR emis- 
sion problem. Our conclusions are summarized in the last 
section. 



2 CONTINUOUS DUST HEATING AND THE 
FOREGROUND EMISSION 

As we mentioned in the introduction, there are two ways to 
tackle the foreground problem. One of them is to assume 
that the optical depth at the darkest spot of the studied 
region is so large that all the emission from this location is 
caused by foreground sources. The other approach is to esti- 
mate the foregroun d emission as a sum of the zodiacal light 
|Kelsall et al.lll9 98) and the interstellar matter contribution 
l|Butler fc Tan 2009). The two approaches may lead to quite 
different inferences. 

Let us take the IRDC 321.73+00 5 (P2) core (IRDC 321 
for short) from IPavlvuchenkov et all (|201ll ) as an example. 
The background intensity at 24 /im in its vicinity is about 
35 MJy/ster, while the intensity at the peak of the millimetre 
emission is less than 30 MJy/ster. The zodiacal light contri- 
bution at 24 /im, as given by the SPOT software, is about 
18 MJy/ster. As this object is relatively nearby (~ 2 kpc, 
IVasvunina et al] 120091 ) , the s moothed Galac t ic for eground 
contribution, estimated as in iButler fc Tanl (2009), is less 
than 10%, so we neglect it here. Subtracting the zodia- 



© 0000 RAS, MNRAS 000, 000-000 



Mid-infrared emission in protostellar cores 3 



cal light only, we get ~ 17MJy/ster for the background 
emission (plus any unaccounted foreground emission) and 
10 MJy/ster for the emission at the m illimetre peak. These 
are in tensities that have been used bv lPavlvuchenkov et al.l 
|201l[ ). Alternatively, if we would assume that all the ex- 
tra emission at the bottom of the shadow is the foreground 
emission, we would end up with zero emission for the core 
and ~ 5 — 7 MJy/ster for the surrounding region. 

To illustrate the si gnificance of the choi c e, we repeat 
here calculations from IPavlvuchenkov et al.l |201lf ) using 
these two approaches to the foreground estimation. Prior to 
giving results we recall the basic equations of the adopted 
RT model with continuous dust heating. The core is assumed 
to be spherically symmetric, with the density distribution 
given by t he same expression as is used for low-mass prestel- 
lar cores |Tafalla et al. II2002F) : 



n(H 2 



no 



(1) 



l+{r/r y 

where no is the central density, ro is the radius of an inner 
plateau, p describes the density fall-off in the envelope. To 
take into account the heating by an embedded protostar (or 
a group of protostars) we assume that there is a black-body 
source in the core centre with the temperature T* and a 
fixed radius of 5 Rq . The inner core radius is set to 50 AU 
(which is the upper limit for the dust sublimation radius), 
and the outer core radius is assumed to be equal to 1 pc. 
The core is illuminated by the diffuse ambient interstellar 
field that is characterized by a certain colour temperature 
Tbg and dilution -Dbg. 

We use the Accelerated Lambda Iteration (ALI) 
method for the radiativ e transfer modelling, which is simi- 
lar to that described in Pavlvuchenkov et alj |2004l ) and in 
Hogerhciide & van der Takl (|20ol T but with modifications 
for thermal radiation (see also Hubenvl 120031 . for general 
principles of the ALI). The mean radiation intensity 



J v = (in)- 1 J I v {n)dQ, 



(2) 



is determined by integrating the radiative transfer equation 
(nV) I v = n„ {S u - Iu) (3) 

along representative directions. Here I v is the specific radia- 
tion intensity, S v is the source function, k v — a v + a v is the 
extinction coefficient, a v is the absorption coefficient, and av 
is the scattering coefficient. The temperature of the medium 
T is found from the equation of radiative equilibrium 



J ct v J v dv — J ct v B v 



(T)du, 



(4) 



where B v is the Planck function. The adopted convergence 
criterion is that the relative difference in temperature be- 
tween subsequent iterations is smaller than 10~ 3 at all radii. 
Although this formalism is multidimensional, in this study 
the RT equation is solved in ID, under the assumption of 
spherical symmetry. 

In the radiative transfer equation, scattering is taken 
into account in the approximation of isotropic coherent scat- 
tering, so that the source function S v takes the form 

a v B v + o v J v 



(5) 



The absorption and scattering coefficients as functions of 
frequency for amorphous silicate grains of a given radius 
are computed using the Mie theory. To calculate the total 
absorption and scattering coefficients for a dust ensemble we 
assume that the size distribution of dust grains is de scribed 
by a power law, f(a) oc a -3,5 (jMathis et al] fl977h . with 
minimum and maximum grain radii of 0.001 and 10 /im. For 
each parameter combination of the core model, we simulate 
the radiative transfer, compute the intensity distributions at 
each considered wavelength, and quantitatively compare the 
computed and observed distributions using the standard \ 2 
criterion. The theoretical distributions are convolved with 
the relevant telescope beams for each wavelength. The search 
for best-fit parame ter values is perform ed with the PIKAIA 
genetic algorithm (|Charbonneaulll995l). More details about 
the fit ting procedure can be found in IPavlvuchenkov et al.l 
i|201ll ). 

Results of the fit are shown in Figure [1] In the left col- 
umn we present data for the case when only DIRBE-based 
zodiacal light estimate is subtracted from the signal (here- 
inafter DZF case). Practically, we subtract zodiacal light 
from 8 um and 24 fim data. Its contribution at 70 (jm and 
1.2 mm is negligible. The right column in Figure [1] corre- 
sponds to the case when the foreground level is assumed to 
be equal to the lowest intensity within the IRDC 321 core 
(hereinafter complete foreground subtraction, CFS case). 

Obviously, in the CFS case we do not have problems ei- 
ther with 24 um emission or with other wavelengths. At the 
same time, some core parameters, that are derived under this 
assumption, are quite different from the results obtained in 
the DZF case. While temperature profiles and central den- 
sities no are almost identical, slopes of the power-law enve- 
lope p axe markedly different (Figure [2j|. In the CFS model 
p = 2.5, which is typical for prestellar objects. In the DZF 
model the density fall-off is much steeper, with p — 3.7. Ap- 
parently, in the DZF case the algorithm tried to reconcile the 
need to have more material to account for emission at longer 
wavelengths and the need to have less material to account 
for low absorption at shorter wavelengths. As a result, the 
core mass is 310 Mq in the DZF model and 750 Mq in the 
CFS model. One may argue that the factor of two difference 
in the core mass is not that significant. However, different 
slopes of the density radial profile in the envelope may lead 
to a greater effect when applied to chemical modelling and 
molecular line RT modelling. 

We should note that in general the CFS model does look 
more attractive. It produces more physical density profile, 
is consistent with the 24 /jm intensity radial distribution. 
Also, analysis of the algorithm convergence shows that the 
solution in the CFS case is well-defined, while in the DZF 
case an alternative solution exists wit h higher central den- 
sity ( ~ 10 9 cm -3 ) and smaller plateau l|Pavlvuchenkov et al.l 
l201lf ) which is nearly as good (by formal x 2 criterion) as the 
low density solution presented here. 

So, after all, to reproduce all the bands, it seems that it 
is only necessary to assume that all the extra emission comes 
from the foreground, i.e. to rely on the CFS approach. But 
how justified is this assumption? Let us consider some gen- 
eral arguments. The distance to the IRDC 321 core that 
we use as an example is about 2kpc. To account for the 
24 jj,m emission, as described above, we would need an in- 
terstellar foreground fraction as high as 70%. However, an 
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Figure 1. Theoretical and observed integrated intensities in the IRDC 321 core computed with different assumptions on the foreground. 
In the left column all the foreground emission is caused by the zodiacal light estimated using the DIRBE model (DZF). In the right 
column it is assumed that the foreground emission intensity is equal to the intensity at the darkest region of the core (CFS). Dots 
represent individual observational pixels. 
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Figure 2. Density and temperature radial profiles for the models 
presented in Figure [l] 



approximate method described in lButler fe Tanl (|2009h gives 
the expected foreground contribution of only about 10% for 
this object. Also, according to th e Galaxy structure map , 
based on the GLIMPSE survey (|Churchwell et al.l 120091) . 
there should not be much intervening material in the di- 
rection of the IRDC 321 core. 

Let us assume that all IRDC cores are identical, that 
is, they all have the same optical depth r and emissivity 
E (at a certain wavelength). Also, let us assume that the 
observed emission G outside of the core projection is the sum 
of the 'true' background emission B (that comes from the 
material behind the core) and the foreground emission F, 
that constitutes the fraction / of the total emission, F — fG 
(Figure [3]). Then, the observed emission contrast 



V = 



G - (Ber T + E + F) 



G 



After some algebra we obtain 
V =(!-/)(! 



E 
G' 



(6) 



(7) 



In Fi gure B] we show the obs ervat ionally inferred r\ for cores 



Vasvunina et all |2009l ) and lButler fc Tanl (|2009h sam- 



from 

pies at 24 /im (after zodiacal light subtraction). As we do 
not intend to make any precise conclusions, values are sim- 
ply read off by eye from the MIPSGAL cutouts. The plot 
apparently shows that r\ grows for G < 30MJy ster -1 and 
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Figure 3. Scheme explaining various symbols entering equations 
@ and {7). 
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Figure 4. Emission contrast 77 at 24 (im for cores from samples 
by Vasyunina et al. (2009) and Butler & Tan (2009) as a function 
of the total IR background intensity. Zodiacal light contribution 
is subtracted. 



then stays constant. If the proper core emission were absent, 
that is, E — 0, we would have r\ independent of G (of course, 
under the assumption that / is independent of G). On the 
other hand, if E is not zero, r\ would be an increasing func- 
tion of G as long as G is not too big to make the second term 
small. This is really observed for G < 30MJy ster -1 . The 
increase of 77 by 0.4 from G = 10 MJy ster -1 to G = 30 MJy 
ster -1 corresponds to E w 5 MJy ster -1 . At larger values 
of G the second term in eq. ([7J) is small, and r\ stays nearly 
constant with some scatter caused by scatter in / (that is, 
in distances and Galactic structure features). 

Solid line in Figure U is drawn using eq. (0) for / = 0.4 
and E = 5 MJy ster -1 (T24 = 00). Thus, observed 77 val- 
ues indicate that there exists a typical foreground frac- 
tio n f ~ 0.4, that is clo se to the value of 0.54 inferred 
bv iPeretto fc Fulled | 20091 ) . This is a statistical value, of 
course, that is not directly applicable to specific objects. For 
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example, in Figure [4] a point corresponding to IRDC 321 
lies above the solid line, which is consistent with / being 
somewhat smaller than 0.4. An analogous plot for an 8 /im 
contrast does not show any clear correlation with G, with 
an average 778 value of 0.6 and significant scatter. This also 
implies / ~ 0.4 (assuming E$ — and rg = 00). 

So, the observed IRDC 24 contrasts are broadly con- 
sistent with the omnipresence of the intrinsic mid-IR emis- 
sion in IRDC cores. In other words, the DZF approach can 
actually be valid, indicating, though, that we miss some im- 
portant mechanism(s) responsible for the intrinsic mid-IR 
emission of the cores. 

It is important to remember that different assumptions 
on the foreground lead to different inferences on the core 
structure, mass, chemical composition, and, by implication, 
on its evolutionary status. If we intend to use the informa- 
tion, deduced from the continuum emission modelling, to 
simulate molecular line profiles, we do need a more substan- 
tiated approach to the foreground problem. In the following 
section we check whether the discrepancy in the 24 fj,m emis- 
sion can be alleviated by taking into account the stochastic 
heating of small dust grains that is known to affect a SED 
in the mid-IR range. 



3 STOCHASTIC DUST HEATING 

To include effects of stochastic dust heating, the RT model 
needs to be upgraded. First, it is no longer possible to as- 
sume that all dust grains have the same temperature irre- 
spective of their size. Second, a small grain is not in ther- 
mal equilibrium with the radiation field, and its temperature 
varies with time. So, in general, dust temperature of a single 
particle is a function of grain size and time. 

To keep the problem tractable, we drop the continuous 
size distribution and assume that dust consists of finite num- 
ber of components, N, so absorption, scattering and emis- 
sion coefficients are given by 



(8) 



Individual absorption and scattering coefficients are given 
by 



i 2,^abs 
i 2 /-.sea 



(9) 
(10) 



where m, cu, Qt,l , an d <9t°i are the number density, radius, 
absorption and scattering efficiency factors for dust grains 
of ith type. 

Since an isolated dust particle is subject to tempera- 
ture fluctuations, we consider a large ensemble of identical 
particles and describe their thermal state by temperature 
probability density distribution P(T). The P(T)dT is the 
fraction of particles with temperatures lying in the interval 
(T, T + dT) . The emission coefficient for each dust compo- 
nent is given by 



3» 



OO 

Ylj P l (T)B v (T)dT, 



(11) 



where P l (T) is the probability density distribution for ith 



component, B V (T) is the Planck function. Distributions 
P l (T) are supposed to be known from the previous itera- 
tion. 

Similarly to the RT model with continuous heating, out- 
lined in the previous section, the method with stochastic 
dust heating is also based on the ALI technique and consists 
of two steps. At the first step, the mean radiation intensity is 
computed at each cell by integrating the RT equation (|3]) in 
ID along representative directions. At the second step, the 
mean intensity J„ is used to update P l (T) for all dust com- 
ponents. Evaluation of P l (T) is based on the Monte Carlo 
method. The temperature evolution of an isolated dust grain 
is simulated and then converted into the probability density 
distribution. A single grain temperature evolution is com- 
puted in three steps: 

• generation of time and frequency sequences of absorbed 
photons; 

• calculation of grain temperature spikes due to absorp- 
tion of photons; 

• calculation of continuous grain cooling due to thermal 
emission between two consecutive absorption events. 

More details about the stochastic heating algorithm can be 
found in Appendix. 

In this study we assume that dust consists of four com- 
ponents: large silicate grains, large graphite grains, small 
graphite grains, and PAH particles. Their parameters are 
listed in TableQ] Absorption and scattering efficiency factors 
Q„ B and Q° ca for si l icate and graphite particles are taken 
from lLaor fc Draine rtl993f). Optical pr operties of PAHs are 
calculated following iDraine fc Lil ((2007), where 30 PAH fea- 
tures (described by the Drude profile) are ta ken into ac- 
count . Heat capacities CV(T) are taken from IDraine fc Lil 
<|200lT) . Simulations of IRDCs are performed using 100 log- 
arithmically spaced radial cells. As we use the Monte Carlo 
method to compute P(T), the resultant distributions are 
quite noisy. So we adopt a much softer approach here, assum- 
ing 5% difference in P(T) between subsequent iterations as 
a stopping criterion. Each model requires about 5 Lambda- 
iterations for convergence. The total computation time is 
about one day per model on an average PC. 

The code was tested against a number of problems. In 
the case of large grains (when stochastic heating is not im- 
portant) results are consistent with temperatures computed 
using the code with continuous dust heating. The stochastic 
heating algorithm was tested separately. In particular, we 
successfully reproduce SEP and P l (T) distri butions shown 
in Figures 2 and 6 from iPopescu et al.l (|201 lh . 



3.1 Stochastic heating by interstellar UV 
radiation 

For the computations below, density distribution and the in- 
ner source p arameters are those from th e best-fit IRDC 321 
model from iPavlvuchenkov et all i|201ll ) with the DIRBE- 
based zodiacal background subtraction. The diffuse inter- 
stellar UV field is represented by a Planck spectrum with 
T bg = 20000 K and dilution of 10" 16 . 

In Figure [5] (left) we present radial dependence of the 
normalized probability density distribution P(T) for PAH 
particles (P(T) for VSGs looks similar). The spread of tem- 
peratures in the core envelope is very broad. While at any 
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Table 1. Dust parameters. 



Component 


Radius, 


Heat capacity pa- 


Density, 


Atomic 


Mass 




cm 


rameter Tj, K 


g cm -3 


mass 


fraction 


Large silicate grains 


3 ■ icr 5 


175 




3.5 


28 


0.70 


Large graphite grains 


2 ■ icr 5 


450 




1.81 


12 


0.15 


Small graphite grains 


3 ■ icr 7 


450 




1.81 


12 


0.05 


PAHs 


7 ■ icr 8 


450 




2.24 


12 


0.10 



* See Appendix for the definition of Td. 



given time most of PAH particles are cold (20 K), their small 
fraction is heated up to 1000 K. Maximum temperatures 
near the protostar also exceed 1000 K. In the rest of the 
core P(T) distribution is relatively narrow, corresponding 
to the equilibrium temperature. 

The SEDs for the model core toward the position, offset 
by 1" from the core centre (to exclude the direct emission 
from the central source), are shown in Figure [5] (right). No 
convolution with a telescope beam is applied in this and 
subsequent SEDs. It is hard to isolate the relative contri- 
bution from various dust components on a single spectrum. 
To show their roles, after simulating the thermal structure 
of the object with all four components simultaneously, we 
computed the emergent spectra separately for each compo- 
nent. More precisely, at the ray-tracing stage we only took 
into account emission, absorption, and scattering by a sin- 
gle dust component. These contributions are shown in right 
panel of Figure [5] with coloured lines. With the solid black 
line we show the combined spectrum computed for all the 
four dust populations. 

Physically, the spectrum can be divided into three in- 
tervals, according to the optical depth shown in Figure[6] At 
mm and sub- mm wavelengths thermal emission from silicate 
and graphite grains dominates the spectrum. The optical 
depth at these frequencies is low (Figure [BJ, and intensities 
depend significantly on details of the thermal structure in 
the interior of the core, in particular, on the presence or 
absence of a protostar. At the IR range (1-24 (jm) thermal 
emission from stochastically heated PAHs is mostly respon- 
sible for the emergent radiation. The total optical depth at 
these wavelengths is relatively high, and the spectrum is not 
sensitive to the presence or absence of the internal heating 
source. At even higher frequencies (A < 1 /im) the emergent 
spectrum is formed by large silicate grains via scattering of 
background radiation. The optical depth here is very high, 
and the spectrum is not affected by the interior of the core. 

In Figure [5] (right) intensities observed toward the cen- 
tre of the IRDC 321 core are indicated by filled circles. Ob- 
viously, we do not solve the mid-IR intensity problem even 
by adding the stochastic heating effects. More accurately, 
if there were only PAHs in the core, combined emission of 
stochastically heated particles in the core envelope and in 
the vicinity of the protostar would provide the needed in- 
tensity at 24 /im. However, in the model with all four dust 
components the emission of PAHs from the core centre is 
completely absorbed by intervening large grains. At shorter 
wavelengths even combined (core+envelope) PAH emission 
is not sufficient to explain observed intensity. 

Apparently, stochastic heating by the standard diffuse 
UV field does not make a noticeable contribution to the mid- 
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Figure 6. Optical depths through the entire core toward the 
position, offset by 1" from the core centre, for various dust com- 
ponents. 



IR spectrum of an IRDC core. The contribution can be made 
more significant by increasing the ambient UV irradiation. 
In Figure [7] we show same plots as in Figure [5] but computed 
with the enhanced ambient UV field, having Tb g = 26000 K 
and dilution of 2 x 10~ 15 . The relative fraction of hot PAH 
particles is much higher in this case, and, as a consequence, 
the IR radiation intensity toward the core projected centre 
(again with 1" offset) also increases, this time matching all 
the observed values. Computed intensities for the millimetre 
emission are higher than observed ones because we do not 
convolve them with the beam size that is quite large at this 
wavelength. 

However, having solved the problem for the central 
spectrum, we have simultaneously created an even more se- 
vere problem for offset positions. This is illustrated in Fig- 
ure [8] where we show distributions of 24 /im intensity for 
the model core with the standard (blue lines) and enhanced 
(red lines) ambient UV fields. With solid lines we show the 
'true' emergent spectrum that includes both the attenuated 
background and the core emission. The proper core contri- 
bution is shown separately with dashed lines. Obviously, for 
the standard UV field stochastic heating makes negligible 
contribution to the overall core emission (blue dashed line), 
and the computed intensity at the projected core centre is 
much smaller than the observed one. As we make the UV 
irradiation stronger, the intensity in the central part does 
grow up to the observed value, but at the same time it gets 
even higher at the core periphery. The radial intensity pro- 
file has a central flat region of about 5" in size and reaches 
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Figure 5. Results of RT simulations for the model with the standard interstellar field. Left: shown with colour is logarithm of the nor- 
malized temperature probability distribution P*(T) for PAH particles (P*(T) = P(T)/ max{P(T)}). Right: spectral energy distribution 
toward the position offset by 1" from the core centre. Spectra are computed separately for large silicates (red lines), large graphite grains 
(thick blue lines), small graphite grains (thin blue lines), and PAH particles (red lines). The emergent spectrum with all dust populations 
taken into account is shown with solid black line. The background interstellar radiation is shown with black dashed line. Filled circles 
indicate observed values for the IRDC 321 core. 
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Figure 7. Results of RT simulations for the model of IRDC 321 with enhanced interstellar radiation field. 



the maximum value at 35" . This ring- like intensity distribu- 
tion is formed by stochastically heated small grains in the 
envelope as schematically shown in Figure [9] The thickness 
of the heated layer is low, so the projection effect plays a 
role in the formation of the intensity radial distribution. In 
other words, column density of heated VSGs and PAHs to- 
ward the centre of the core is lower than toward the core 
periphery, which results in the intensity difference. The in- 
tensity distribution at 70 /im over the core surface is similar 
to the one for 24 except for the strong emission peak 
toward the location of the protostar. This peak appears due 
to lower optical depth at 70 /im. 

Such ring-like structures, apparently, are not ubiquitous 



in IRDCs, so we must admit that stochastic heating do not 
seem to be a viable explanation for the alleged excess mid-IR 
emission in IRDC cores. But the problem would be solved 
if we can find some other dust heating mechanism which 
would work not only in the core envelope but in its whole 
volume. Among various non-radiative dust heating mecha- 
nisms, there is one that is often included in chemical models 
of prestellar cores. Specifically, these models include man- 
tle evaporation due to stochastic hea ting of dust grains by 
cosmic ray particles ( Leger et al 1 19851 ). In the following sub- 
section, we check if same dust heating that is responsible for 
mantle destruction can change appreciably the core SED. 
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Figure 8. Intensity distributions at 2&fim over the model core 
for cases with and without the background radiation as well as 
with and without UV stochastic heating. Observed intensity for 
the IRDC 321 core is schematically shown with gray band (DZF 
case) . 
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Figure 9. Scheme explaining the formation of the ring-like inten- 
sity profile in mid-IR for the core heated by interstellar radiation. 

3.2 Stochastic heating by cosmic rays 

The interaction of cosmic rays (CR) with interstellar 
medium is a complex process. The collision of highly en- 
ergetic CR particles with interstellar molecules results in 
formation of secondary particles (pions, mesons, positrons, 
gamma-rays, etc.) that also interact with interstellar gas. 
Both primary and secondary non-thermal particles dissoci- 
ate and ionize interstellar molecules as well as collide with 
dust grains and contribute to their stochastic heating. Ide- 
ally, to calculate stochastic heating of dust grains induced 
by cosmic rays one should know 1) the flux/energy distri- 
bution of all non-thermal particles, and 2) the fraction of 
the particle kinetic energy that goes into the grain thermal 
energy upon the collision. The proper modelling of these 



processes is a challenging problem that is further restricted 
by uncertainties of involved physical data. 

Here we use a simplified approach. The most significant 
contribution to stochastic heating of dust grains is assumed 
to come from free electrons formed due to propagation of 
cosmic rays. The mean energy of free electrons produced 
by co smic ray ionization is 20-35 eV (Glassgol d fc LangerJ 
1 97.'! ) . We assume that this energy entirely goes to the ther- 
mal energy of a dust grain in a single collision. We adopt the 
standard ionization rate 10 _17 s _1 and assume that these 
electrons do not interact with gas and only collide with dust 
grains. Since the last assumption is unrealistic (the total ef- 
fective cross-section of molecules is much higher than the 
total effective cross-section of dust grains), we significantly 
overestimate the effect. 

In Figure [TO] we present P(T) distribution for PAH par- 
ticles in the model with CR stochastic heating. To isolate 
the CR effect, we treat heating by external UV radiation 
in the integral non-stochastic way which provides the min- 
imum temperature for grains. While PAHs are sometimes 
heated up to 1000 K, the probability to find a particle with 
such a temperature is low. In other words, heating events 
are very rare, and PAHs spend most of the time at the low- 
est temperature. The contribution of grains, stochastically 
heated by CR, to the SED is even smaller than in the case 
of standard UV (right panel of Figure I10[) . The computed 
24 /im intensity is only a few hundredth MJy/ster, while 
the observed intensity is about twenty MJy/ster. Given the 
significant overestimation of the collision rate between elec- 
trons and grains in our model, we conclude that stochastic 
heating of dust grains due to CR cannot solve the problem 
of high mid-IR intensity toward IRDC cores. 



4 DISCUSSION 

The analysis of IRDC mid-IR properties indicates that these 
clouds may possess intrinsic emission in this range. On a 
larger scale, it is typically assumed that emission in the 
24 /im Spitzer band is produced by stochastically heated 
small grains. Our study shows, however, that on a smaller 
scale stochastic dust heating by UV photons and CR parti- 
cles cannot solve the p roblem of excess 24 /im emission, that 
iPavlvuchenkov et all <|201ll ) have inferred for some IRDC 
cores. 

For a standard interstellar UV field the contribution 
of stochastically heated grains to the core mid-IR emission 
is negligible. If we try to overcome this problem by enhanc- 
ing the ambient UV field, the intensity distribution becomes 
ring-like. The reason behind this is quite simple and qual- 
ita tively similar to th e "limb brightening" effect considered 
bv lLeung et~aH (jl989h for longer wavelengths (140-300 (jm) 
which are not affected by stochastic heating of small grains. 
The UV irradiation of the core creates the overheated dust 
layer on its surface. Geometrically, this layer is narrow due 
to strong UV absorption. On the other hand, it is nearly 
transparent in the mid-IR range so that the hot surface of 
the dense core should be observed as a bright rim. Such 
structures are not seen in any of the cores we have checked 
in this study. While 24 pm rings are actually observed in 
other sources, most of them are classified as planetary neb- 
ulae and circumstellar shells around massive (proto) stars 
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Figure 10. Same as in Figure [5] but for the model with CR stochastic heating. 



ijMizuno et al.ll2010l : IWachter et"aDl201dl ^ . We conclude that 
the stochastic heating of small particles, at least, in its typ- 
ical implementation, does not seem to be a viable candidate 
for the source of mid-IR emission in IRDC cores. Other 
known potential stochastic heating factor, cosmic rays, is 
also not effective, even if somewhat extreme parameters are 
adopted. 

Our inability to fit mid-IR observations may at least 
in part be related to the deficiencies of the model. For ex- 
ample, we assume spherical symmetry for our cores, with 
smooth density distributions, resulting in large UV optical 
depths in the core outer parts. The problem with 24 /im 
emission would probably be alleviated if UV radiation is 
able to penetrate deeper into the core due to its irregular 
structure. In this case instead of the narrow hot layer we 
would have a more extended hot envelope. Also, we consider 
only UV and CR heating, but a nother bulk emission gen - 
eration mechanisms are possible. iDulev fc Williams! |201ll ) 
suggested that surface chemical reactions can be an energy 
source that excites infrared emission bands even in grains of 
'classical' sizes. Some chemically induced energy deposition 
in dust grains is also implied by t he non-therma l mant le des- 
orption mechanism sugg ested by iGarrod et all (|2007h . The 
presence of some unaccounted energy source indirectly fol- 
lo ws from the low sticking pro bability (of about 0.3) found 
bv lPavlvuchenkov et al.l l|2006l ) in their detailed study of the 
CB17 core. They suggested that this low value may be an 
indication of some desorption mechanism which is present in 
these objects and is different from thermal desorption and 
cosmic ray induced desorption. Such a desorption mecha- 
nism would also mean an energy deposition into dust grains. 

If the core irregular structure is related to turbulence, 
then turbulent dissipation may provide some heating. How- 
ever, it seems likely that in the absence of strong shocks this 
mechanism would rather somewhat raise the mean temper- 
ature of grains, enhancing the core emission in the far-IR 
and submillimetre bands, but not in the near- and mid-IR. 

Yet another factor that is capable to affect the mid-IR 
core emission is scattering. In this study, we assume that 



scattering is isotropic. If our assumption is relaxed in favour 
of forward scattering we would have some additional con- 
tribution from the core central hot region. However, the in- 
ferred core structure implies that it is opaque at 24 /im, so 
even with forward scattering the contribution from the inner 
region should not be important. 

Scattering of the ambient infrared emission can be 
significant if we assume some grain evolution in the 
cores. Extended mid-IR emission in the 3.6 and 4.5 /tm 
bands, that was terme d MI R cloudshine and co reshine by 
ISteinacker et al.1 (|201Ch and iPagani et alj < |2010h . observed 
in some molecular clouds and cores, was interpreted as a 
result of scattering of the ambient infrared light by large 
grains. In principle, non-zero E due to scattering at 24 fim 
is also possible. In this study we consider only four rep- 
resentative grain types and do not vary their parameters 
along the core radius. If we would include very large grains 
(a > 10 /im) in the model, scattering at 24 fj,m would be 
greater, making mid-IR shadows less prominent. However, 
these large grains would need to be present in sufficient 
amount such that scattering dominated over the absorption 
by small grains. It is also possible that scatte ring on large 
porou s dust aggregates, like those considered bv lOrmel et al.l 
may contribute into observed mid-IR flux. 

Obviously, by varying dust optical properties and size 
distribution along the core radius, we may achieve a better 
fit. However, these parameters cannot be varied in a uncon- 
strained way. Grains in the inters tellar medium can grow by 
accre tion of refractory elements jVoshchinnikov k, Hennina 
I2010I) or volatile species but the potential for this growth 
is limited by the abundance of heavy elements. So, the pre- 
ferred way for grains to grow to large sizes is coagulation. In 
other words, the number of large grains can only increase at 
the expense of small grains and PAH particles. The two-fold 
outcome of this process is that scattering becomes more ef- 
fective, but at the same time contribution of stochastic heat- 
ing diminishes. To estimate the physically motivated balance 
between these two processes, one needs to include a grain 
growth in the model. In general, the search for dust pa- 
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rameters is a complex problem which has to be solved in 
a self-consistent way, using observations at as many wave- 
lengths as possible. Such a modelling should be a subject of 
a dedicated study. 



5 CONCLUSIONS 

Recently, IPavlvuchenkov et al.l |201ll ) found that 24 fim 
emission in mid-IR shadows of typical IRDC cores does not 
fit into the combined emission and absorption picture, be- 
ing too bright for the density and temperature distributions 
that are needed to explain intensity maps of the cores in 
other IR and millimetre bands. This problem arises if the 
only foreground that is subtracted from the signal is the zo- 
diacal light contribution estimated from the DIRBE data. 
The natural way out of this is to assume that the remaining 
'extra' emission comes from some other foreground source, 
like intervening interstellar dust. However, analysis of obser- 
vational data shows that the typical foreground contribution 
is about 40% or less for nearby objects. All the remaining 
emission should originate in the object itself. 

We checked if stochastic heating of small grains can 
produce this extra 24 /im emission. However, it only allows 
to reproduce the central mid-IR intensity for the studied 
core if we assume that the core is illuminated by the diffuse 
UV light enhanced by a factor of 70 relative to the aver- 
age Galactic value. In this case, intensity at the core edges 
also grows due to projection effects so that the entire core 
acquires a ring-like appearance which does not seem to be 
observed. Stochastic heating by cosmic ray particles does 
not change the SED significantly due to the low energy den- 
sity input of these particles. 

Overall, we conclude that the origin of at least some 
mid-IR emission in IRDC dark cores is unclear and requires 
a special study. A great care must be taken when someone 
uses infrared observations to deduce IRDC properties and 
to compare them with results of numerical modelling. The 
assumption that IRDCs are completely dark in the infrared 
(so that all the emission is the foreground emission) does 
allow to reproduce the observational data, but in this case 
all the information on the proper core emission is lost which 
may lead to wrong inferences. 
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APPENDIX A: DETAILS OF STOCHASTIC 
HEATING ALGORITHM 

Here we provide details of the method that is used to cal- 
culate the temperature evolution of an isolated dust grain 
exposed to the radiation field with the mean intensity J„. 
Emission (radiative cooling) of the dust grain is described 
as a continuous process, while absorption (radiative heating) 
is represented either as a discrete process or as a continu- 
ous process, depending on the photon energy. The concept 
of contin uous cooling is a r easonable approximation to the 
problem (|Draine fc Lill200l[) . 

First, we split the spectrum of the external radiation 
field into the low-energy and high-energy intervals sepa- 
rated by a critical frequency v c , defined by the relation 
hv c — 0.01-Bth, where E t h is the mean thermal energy of 
the dust grain calculated under the assumption of continu- 
ous heating and cooling. We use stochastic treatment only 
for photons from the high-energy interval of the spectrum 
since low-energy photons does not produce significant tem- 
perature fluctuations of the grain. 

The next step is to evaluate the time and frequency 
sequence of the absorbed photons within the interval 
(^c,^ max ), where f ma x = 10 16 Hz is the adopted maximum 
frequency of photons. We consider M absorption events 
and generate the frequency sequence \y\, vm} using the 
Monte Carlo simulation for the re-normalized probability 
density distribution of the absorbed photons 



Qabs t 

hv 



(Al) 



where Q„ bB is the absorption efficiency factor for a given 
grain type. The corresponding time sequence {ti, tjvf } of 
the absorption events is simulated using the Poisson statis- 
tics 



f{t) = Xexp(-Xt), 



(A2) 



where f(t) is the probability density distribution for the 
time interval t between successive events, A = M/to is the 
mean number of events per unit time, to is the total length 



of the sequence. The total time to can be determined from 
the relation 



f abs 1 
io I Ql B J v dv = — ^2 hv 



(A3) 



which represents the energy absorbed by the grain. The set 
of obtained sequences {i/i} and {ti} is the discrete represen- 
tation of the energy deposited into the grain during time to 
by photons from (i/ c ,i/ max ). 

Next we calculate the temperature jump of the grain 
due to the photon absorption. Let us suppose that the tem- 
perature of the grain just before the absorption is To. The 
temperature Ti just after the absorption of a photon with 
energy hv is given by 



U(Ti) - C/(T ) = hv, 



(A4) 



where U(T) is the thermal energy of the grain with temper- 
ature T. The relation between thermal energy and temper- 
ature is non- linear since the heat capacity of the grain, Cv, 
is a function of temperature. We approxima t ed hea t capac- 
ities, presented in Figure 2 of iDraine fc "LI (|200lT ). by the 
simple phenomenological law 

" , (A5) 



Nk 



1 + 



where N is the number of atoms in the grain, k is Boltzmann 
constant, the heat capacity parameter Td is 175 K for silicate 
and 450 K for graphite. The corresponding thermal energy 
as a function of temperature 



U 
Nk 



3 [T — Td arctan 



(A6) 



Equation (|A6f) is substituted into equation (|A4|) . which is 
then solved for Ti using the bisection method. 

In order to calculate the temperature evolution of the 
dust grain between two subsequ ent absorption events we 
solve the following equation (e.g.. iKruegell [20031 ): 



dt 




I'm ' 

/ 



J v dv- / Ql B v {T)dv 



(A7) 

The first term on the right hand side is the energy absorbed 
from the radiation field in the low-energy frequency inter- 
val per unit time, while the second term is the cooling rate 
due to continuous radiation. We solve this equation using 
an implicit Euler method where the corresponding finite- 
difference equation is solved by the bisection method. Since 
the solution of this equation describes a fast temperature 
decay right after the absorption followed by the slow tem- 
perature evolution, we use adaptive time step control which 
allows having about 30 time-steps to evaluate the tempera- 
ture evolution between absorptions. An example of the tem- 
perature history for a small graphite grain is shown in Fig- 
ure [Si] 

The last step of the algorithm is to convert the history 
T(t) into the temperature probability density distribution 
P(T) assuming the ergodic hypothesis. We split the temper- 
ature into the number of intervals and calculate the relative 
time that the grain spends in each interval. An example of 
P(T) for different grain sizes is shown in Figure [ 
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Figure Al. The temperature evolution of a 10 — 6 cm graphite 
grain exposed to the Planck radiation with T^ s = 20000 K and 
dilution of 10" 16 . 
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Figure A2. The temperature probability density distribution for 
graphite grains of different radii exposed to the Planck radiation 
with T bg = 20000 K and dilution of 10~ 16 . 
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